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In the present work we investigate the adequacy of broken-symmetry unrestricted density func- 
tional theory (DFT) for constructing the potential energy curve of nickel dimer and nickel hydride, 
as a model for larger bare and hydrogenated nickel cluster calculations. We use three hybrid func- 
tionals: the popular B3LYP, Becke's newest optimized functional Becke98, and the simple FSLYP 
functional (50 % Hartree-Fock and 50 % Slater exchange and LYP gradient-corrected correlation 
functional) with two basis sets: all-electron (AE) Wachters+f basis set and Stuttgart RSC effective 
core potential (ECP) and basis set. We find that, overall, the best agreement with experiment, 
comparable to that of the high-level CASPT2, is obtained with B3LYP/AE, closely followed by 
Becke98/AE and Becke98/ECP. FSLYP/AE and B3LYP/ECP give slightly worse agreement with 
experiment, and FSLYP/ECP is the only method among the ones we studied that gives an unacept- 
ably large error, underestimating the dissociation energy of M2 by 28 %, and being in the largest 
disagreement with the experiment and the other theoretical predictions. We also find that for M2, 
the spin-projection for the broken-symmetry unrestricted singlet states changes the ordering of the 
states, but the splittings are less than 10 meV. All our calculations predict a 55-hole ground state for 
N12 and S- hole ground state for NiH. Upon spin-projection of the singlet state of M2, almost all of 
our calculations: Becke98 and FSLYP both AE and ECP and B3LYP/AE predict 1 {d^ 2 _ y2 d^ 2 ^ y2 ) 

or 1 (d xy d^ y ) ground state, which is a mixture of J E+ and 1 T g . B3LYP/ECP predicts a 3 (d^ 2 _ y2 d^ y ) 
(mixture of 3 E~ and 3 r„) ground state virtually degenerate with the 1 (d^ 2 - y zd^2 _ y2 ) f 1 (d^ y d^ y ) 
state. The doublet 5-hole ground state of NiH predicted by all our calculations is in agreement 
with the experimentally predicted 2 A ground state. For M2, all our results are consistent with the 
experimentally predicted ground state of 0^J" (a mixture of 1 Ej" and 3 Eg ) or 0„ (a mixture of 1 E^" 
and 3 E+). 



I. INTRODUCTION 

During the last decades clusters have been extensively 
studied because of their potential applications, their the- 
oretical value in understanding the transition from iso- 
lated atomic systems to condensed matter 0, Q and their 
relevance to the study of surface processes and heteroge- 
neous catalysis The rapid development of exper- 
imental techniques in recent years has made it possible 
both to obtain size-controlled transition metal clusters 
and to study their reactivity against chemisorption pro- 
cesses USE 13 

Methods for studying properties and behavior of clus- 
ters have been developed, and a review on computa- 
tional studies of clusters has been written by Freeman 
and Doll There have been many studies on nickel 

clusters using various methods of exploring the potential 
energy surfaces (PES). The construction of such poten- 



* Electronic address: lcvdiaconu@brown.edul 

T Present address: Department of Chemistry and Center for 
Biomolecular Simulation, Columbia University, New York, New 
York 10027. 



tial surfaces is a major problem, especially for transition 
metal clusters. Many methods have been used to con- 
struct PES's for nickel clusters, ranging from empirical 
- Finis-Sinclair type [H HI], semi-empirical — tight 
binding 0, Q] and extended Hiickel [l^, to ab initio 
or mixed empirical-a& initio [16| approaches. Recently, 
there have been studies of hydrogen atoms on Cu sur- 
faces 17] within the density functional framework. It 
has been found that semiempirical methods are insuffi- 
cient for accurate description of such systems, and first 
principle quantum-mechanical methods are needed to ob- 
tain a proper description of the hydrogen binding site. 

Our long-term goal is to explore the structure and dy- 
namics of clusters, including nickel and nickel hydride 
systems. The combination of the physical complexity 
and the computational demands of these systems neces- 
sitate that the microscopic force laws that are utilized in 
such simulations be both efficient and reliable. 

Among the correlated electronic structure methods the 
best candidate is clearly density functional theory (DFT) 
because of its ability of reaching high accuracy — similar 
to coupled-cluster CCSD(T) method for second-row ele- 
ments — when hybrid exchange-correlation functionals 
are used ^3- Moreover, DFT (using hybrid function- 
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als) is computationally not much more expensive than 
Hartree-Fock. 

While there have been a number of DFT calculations 
reported on small nickel clusters [H S3, I2J |H SHH 
l25l l2fi| , some results appear to be inconsistent both with 
respect to available experimental data and/or with re- 
spect to other theoretical predictions. 

The works of Yanagisawa et al. j2^| and Barden et 
al. [2^1 on the performance of DFT on the first transi- 
tion metal series have shown that non-hybrid functionals 
(BLYP SI H3, BP86 |H |H H2, BOP (2lj3i and 
PW91 _0ia nd hybrid functionals (B3LYP |2fll2a, BH- 
LYP |3(l |3J|) give an overall similar description for 3d 
transition metal dimers, with the non-hybrid ones giving 
better bond lengths and the hybrid ones better dissocia- 
tion energies. However, while Yanagisawa et al. |27| ob- 
tain good agreement with experiment for nickel dimer for 
all studied exchange-correlation functionals (they only 
calculated the triplet states), Barden et al. |28j obtained 
a negative dissociation energy for their calculated singlet 
ground state with the B3LYP functional (and negative 
or very close to zero for all hybrid exchange-correlation 
functionals). This prompted us to use symmetry break- 
ing in unrestricted DFT for describing the lowest singlet 
state of nickel dimer. With larger cluster calculations in 
mind, we also used broken symmetry unrestricted DFT 
to better describe bond breaking in all states of nickel 
dimer and nickel hydride. 

It has been argued that broken-symmetry unrestricted 
calculations (Hartree-Fock and DFT with hybrid func- 
tionals) are useful for describing systems with weakly 
coupled electron pairs EE E3, EE I^S E3, EH E2- Ni 2 
is definitely such previously observed by Basch 

et al. 19]. As argued by Cremer [lq. the combination 
of hybrid exchange-correlation functional with symmetry 
breaking leads to a better description of systems in which 
static correlation is present than does the restricted DFT 
formalism. Finally, we believe that the formalism used 
to describe any system is solely dictated by the objec- 
tive of the calculation. For a variational approach, and 
DFT can be regarded as such — aside from the exchange- 
correlation functional, — the more flexible is the form of 
the trial function (density), the lower is the obtained en- 
ergy. Since our interest is mainly in the energetics of 
nickel clusters, the best choice for us seems to be the un- 
restricted broken symmetry DFT approach with hybrid 
functionals. 

In the present work we study the nickel dimer and 
nickel hydride using broken symmetry unrestricted DFT 
with hybrid exchange-correlation functionals — mainly 
the popular B3LYP j33, |3j| - - as model systems for 
larger bare and hydrogenated nickel clusters in an at- 
tempt to establish what might comprise a minimally re- 
liable method for more extensive nickel cluster calcula- 
tions. 

The outline of the reminder of the present paper is as 
follows: SectionlTTlcontains discussions the methods used, 
Section ITTT1 presents the results of the calculations, and, 



where possible, comparisons with previous reports. Sec- 
tion I1VI concludes with suggestions for further research 
based on the present findings. 

II. COMPUTATIONAL DETAILS 

The DFT calculations reported in this paper are car- 
ried out with NWChem [43 computational chemistry 
package, using the unrestricted Kohn-Sham approach, al- 
lowing for symmetry breaking, and using a finite orbital 
(spherical Gaussian) basis set expansion and charge den- 
sity fitting. 

Hartree-Fock (HF) and second order M0ller-Plesset 
(MP2) calculations are performed for comparison for 
the states of the Ni atom and are done in unrestricted 
form.fr^ 

Throughout the paper we will use the notation 
M (h h B ) for the states of nickel dimer, where M is the 
multiplicity, h A and h B are the unoccupied (hole) or- 
bitals in the 3d shell on the two Ni atoms, denoted A and 
B. The broken symmetry singlet states (with S z = and 
(S 2 ) = 1) are denoted by 1 ^{h A h B ). 

In general, an unrestricted Slater or Kohn-Sham deter- 
minant is not an eigenfunction of the total spin operator 
S 2 , and the results can only be characterized by the num- 
ber of a and (3 electrons. However, following common 
usage, we refer to the states that differ in the number 
of a and j3 electrons by as singlets, by 1 as doublets, 
and so on. We explicitly identify pure spin states where 
relevant. 



A. Exchange-Correlation Functionals 

We used three hybrid exchange-correlation function- 
als: the very popular B3LYP — composed of the B3, 
Becke's three-parameter hybrid exchange functional [35| 
and LYP |3(j correlation functional — is the first choice 
because it is well known and extensively characterized. 
Becke's newest optimized functional, Becke98 Q] is also 
used, since it is supposed to be, in a certain sense, the 
best obtainable exchange-correlation functional within 
the gradient-corrected framework. The hybrid composed 
of half Slater exchange 0] , half Hartree-Fock exchange 
and LYP 30] correlation, named here FSLYP is also 
used for comparison, as it is the simplest theoretically- 
justifiable hybrid method and is reported to perform 
rather well [SliJ]. 

B. Basis sets 

All calculations are performed with spherical basis 
sets. As all-electron (AE) basis sets, Wachters+f ba- 
sis set mm m S3, a ^^p^^i^p^f) c ° n - 

traction is used for nickel and 6-311++G(2d,2p), a 
[6s2p]/(4s2p) contraction for hydrogen. 
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Effective-core potentials (ECP) are also explored, since 
they greatly reduce computational cost. Stuttgart RSC 
ECP effective core potentials basis set [^3, are use d 
for nickel, as they provide a similar quality of valence 
basis functions as Wachters+f. 

Ahlrichs Coulomb Fitting jSJ basis is used as a 
charge density (CD) fitting basis only for the all-electron 
calculations, as it significantly reduces computing time, 
especially for larger systems. When not specified other- 
wise, all reported all-electron results are obtained using 
charge density fitting. 

We did not use charge density fitting with ECP be- 
cause of the large errors that resulted when we tried the 
use of Ahlrichs Coulomb Fitting basis in combination 
with Stuttgart RSC ECP. For example, for B3LYP func- 
tional, CD fitting error is as much as 0.3 eV for both the 
interconfigurational energies of Ni atom and the binding 
energy of Ni 2 . Please refer to Appendix [S] for discussion 
of the accuracy of charge density fitting. 



C. Numerical integration and convergence 

The numerical integration necessary for the evalua- 
tion of the exchange-correlation energy implemented in 
NWChem uses an Euler-MacLaurin scheme for the radial 
components (with a modified Mura-Knowles transforma- 
tion) and a Lebedev scheme for the angular components. 
We use three levels of accuracy for the numerical inte- 
gration that are used in our DFT calculations, labeled 
by the corresponding keywords from NWChem (medium, 
fine and xf ine). 

The reported atomic calculations are those obtained 
with the xf ine grid. For geometry optimization and vi- 
brational frequency calculations we use the fine grid. 
And for the potential energy curve (PEC) scans we used 
the medium grid. The maximum number of iterations is 
set to 100 in all calculations. 

Please refer to Appendix [B] f° r details on numerical 
integration and convergence criteria. 



D. Initial guess 

For all DFT methods we first performed a calcula- 
tion for Ni atom using fractional occupation numbers 
(FONs) [54|], as implemented in NWChem. We use an ex- 
ponent of 0.01 Hartree for the Gaussian broadening func- 
tion. We then use the molecular orbitals from the FONs 
calculation, after proper reordering, as initial guess for 
computing the 3 F and 3 D states of Ni atom. We use 
t% e g configuration in the Oh symmetry group for the 3 F 
state. In order to obtain the lowest energy possible for 
the 3 D state, we scan all hole positions: d z 2, d x 2_ y 2 and 
d xy using D^h symmetry group, and d xz and d yz using 
Dih symmetry group, enforcing the position of the hole 
with a maximum overlap condition. 



For Ni 2 and NiH, we use a broken-symmetry initial 
guess of the form: 3d 9 4s 1 + J, J, 3d 9 4s 1 for singlet Ni 2 , 
3d 9 4s 1 tt + IT 3d 9 4s 1 for triplet Ni 2 and Ni3d 9 4s 1 |T 
+ I Is 1 H for NiH. As initial guess molecular orbitals 
we use those from the Ni atom calculations, sweeping 
through all unique positions of the holes in the 3c? orbitals 
of Ni atom(s), and enforcing the position of the hole(s) 
with a maximum overlap condition. 



E. Geometry optimization 

Geometry optimizations are performed using the 
DRIVER module of NWChem using NWChem's default 
convergence criteria (in atomic units): 4.5 • 10~ 4 maxi- 
mum and 3.0 • 10~ 4 root mean square gradient, 1.8 • 10~ 3 
maximum and 1.2 • 10 -3 root mean square of the carte- 
sian step. These convergence criteria give a maximum 
error in equilibrium bond length of less than w 10 -3 A 
for Ni 2 and less than w 5 • 10" 4 A for NiH. The available 
precision is set to 5 • 10 -7 Hartree for the fine grid and 
5 • 10™ 8 Hartree for the xfine grid. 



F. Vibrational frequencies 

Harmonic vibrational frequencies are calculated using 
NWChem's VIB module with the default options. Since 
analytical Hessian for open shell systems is not available 
for the exchange-correlation functionals used, the Hessian 
is computed by finite differences with A = 0.01 Bohr, 
which gives an estimated error for the vibrational fre- 
quencies of s» 0.5 cm -1 (ss 0.25%) for Ni 2 and « 2 cm -1 
(« 0.1%) for NiH. 

G. Spin and Symmetry Projection 

In general, an open-shell Slater or Kohn-Sham determi- 
nant is not an eigenfunction of the total spin operator S 2 . 
However, spin-adapted configurations can be obtained as 
combinations of (a small number of) restricted determi- 
nants |55l l56j| . Unrestricted determinants are not eigen- 
function of the total spin operator S 2 , either, and they 
cannot be spin-adapted by combining a small number 
of unrestricted determinants |55|. However, for antifer- 
romagnetic coupling of two weakly interacting identical 
high spin monomers, Noodleman 3g| derived an approx- 
imate spin projection scheme that is correct to the first 
order in the overlap integrals. Ni 2 can be well approxi- 
mated by such a model. 

As previously observed by Basch et at [lj|, the elec- 
tronic structure of nickel clusters corresponds roughly to 
a model in which the 3c? electrons can be viewed as weakly 
interacting localized 3c? 9 units bound together primarily 
by 4s electrons. If the 4s electrons are paired in a a 
bond, then Ni 2 has two possible spin states: singlet and 
triplet. However, the open-shell singlet state can not 
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be represented by a single determinant, and the broken- 
symmetry single determinant if 3 obtained by putting 
one of the open-shell electrons in a spin a d orbital on 
one of the Ni atoms and the other electron in a spin (3 d 
orbital on the other Ni atom is not pure singlet, but an 
equal mixture of singlet and triplet (using \S,S Z ) nota- 
tion for the spin states): 

* B = ^|0,0> + i=|l,0) 

with the expectation value of the total spin 
b\S 2 \^ b) = 1- In agreement with this model, 
for the broken symmetry calculations of the .!5 Z = 
state of the Ni dimer the expectation value of the 
total spin (S 2 ) is close to the exact value of 1 for 
the broken-symmetry mixed state, and for the triplet 
(S z = 1) state, (S 2 ) is close to the exact value of 2 (in 
both cases, the relative absolute differences between 
the computed and the exact values are less than 2%). 
Mullikcn population analysis also supports the weakly 
interacting 3d 9 units model. For the triplet nickel dimer 
there is a Mullikcn spin population of 1.00 on each Ni 
atom, and for the broken symmetry singlet there is a 
Mulliken spin population of 1.1 on one of the Ni atoms 
and —1.1 on the other. 

Using the approximate projection method of Noodle- 
man [33, the energy of the pure singlet state, -E'(O) can 
be obtained from the energy of the unrestricted broken- 
symmetry singlet, Eb , and the energy of the triplet, E(l) : 

E(0) = 2E B -E(1). (1) 

The same result can be also obtained by the spin projec- 
tion technique (see, e.g., Refs. 1571 and lH^) . 

M2 belongs to D^h point symmetry group, and the 
irreducible representations (irreps.) are good quantum 
numbers for the molecular states. We combine the spin 
projection with symmetry projection to extract the max- 
imum information possible from the single-determinant 
Kohn-Sham DFT calculations. From simple group- 
theoretical considerations one can find that the pure spin 
and symmetry states of N12 that arise from d$ orbitals, 
which arc found to give the lowest energy states for all 
calculations, are: 1 E+, 'Tg, 1 E~, 3 E~, 3 E+ and 3 r„. 
Within the model of two weakly interacting 3d 9 units, 
for the purpose of projection we consider only the active 
electrons and the active orbitals on each center, namely 
d A 2 _ y2 , d A y , d B 2 _ y2 andd*. 

The projection has been carried out using the pro- 
jection operators technique in D&h, which the smallest 
subgroup of Dooh in which all irreps. arising from the 
{d A ) 1 {d B ) 1 configuration can be completely correlated, 
and the following equations relating the energies of the 
pure spin and symmetry states listed above to the ener- 
gies of the computed triplet and projected singlet states 
are obtained: 



(\d A _ y2 d B 2 _y 2 )) = 


1 
2 


[E 


(% + ) 


+ E(%)] 


(2a) 


(^(d A 2 _ y2 d B 2 _ y2 )j = 


1 
2 


[E 




+ E( 3 T U )] 


(2b) 


^C(«)) = 


1 

2 


[E 


(% + ) 


+ E( 1 r g )] 


(2c) 


E( 3 (d A yd B y)) = 


1 

2 


[E 




+ E( 3 T U )} 


(2d) 


E(\d A 2 _ y2 d B y )) = 


1 

2 


[E 


(X) 


+ E( 1 r g )] 


(2e) 


E( 3 {d A _ y2 d B y j) = 


1 
2 


[E 


(%-) 


+ e( 3 t u )}. 


(2f) 



These equations contain the maximal information that 
can be obtained from single-determinant calculations. 

From equations © we can derive the (partially) 
symmetry-adapted equivalent of Eq. Q: 

E (\h A h B )) = 2E {^(h A h B )) - E ( 3 (h A h B )) (3) 

where (h A h B ) represents each of (d A 2 _ y2 d B 2 _ y2 ), 
{d A 2 _ y2 d B y ), and (d A y d B y ). The spin projection has to 
be done separately for each of the combinations of holes 
(d\ 2 d B 2 2 ) and (d A 2 2 d B ). 

\ x z —y z x z — y z l V x z —y z xy / 

Since the equations for the states M (d A y d B y ) have a 
similar form to those for the M (e?A, 2 d B 2 . 2 ) states, 

x y x y 

M {d A 2 _ y2 d B 2 _ y2 ) and M (d A y d B y ) states should have the 
same energy (M can be 1, 3 or (1,3)). We calculate the 
1 ' 3 (d A y d B y ) and 3 (d A y d B y ) states for consistency check. 

Since the bond lengths for the pure spin states are 
different from each other and from the mixed state, we 
use a harmonic approximation of the potential around 
equilibrium bond length for each state: 

E(d) - -D e + ^iu 2 e (d-d e ) 2 

and solve the resulting equations for d e (equilibrium bond 
lenth), D e (dissociation energy) and uj e (vibrational fre- 
quency) for the projected state (here [i denotes the re- 
duced mass of the molecule). 

III. RESULTS AND DISCUSSIONS 
A. Nickel Atom 

The ground state of the nickel atom is 3 F 4 (3d s 4:S 2 ) |5^, 
IfSlj . However, since our calculations do not include spin- 
orbit coupling, we use weighted averages over the J com- 
ponents of the experimental data for comparison, which 
makes 3 D(3d 9 4s 1 ) the ground state, with 3 F(3d 8 4s 2 ) 
state only 0.03 eV higher, and 1 S(3d 10 ) state 1.74 eV 
above the ground state. 

As first estimated by Martin and Hay [(5(j and con- 
firmed by full relativistic calculations done by Jeng and 
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TABLE I: Energies of atomic states of Ni. Values are in eV, relative to the ground state. 











UHF 


MP 2 


FSLYP 


B3LYP 


Becke98 


State 


Exp. ° 


RC b 


Exp.-RC c 


AE 


AE 


AE 


ECP 


AE 


ECP 


AE ECP 


3 £)(3d 9 4s 1 ) 
3 F(3d 8 4s 2 ) 
^d 10 ) 




0.03 
1.74 


-0.36 
0.21 




0.39 
1.53 


1.44 


5.81 


0.27 
1.41 



0.12 


2.62 


0.32 


3.02 




0.36 
1.90 


0.01 


2.21 


0.20 
0.29 
1.78 2.37 



"Weighted averages over the J components of the experimental valu es l59l . 
b Martin and Hay estimations of relativistic corrections from Ref . l6Ct 
c Experimental values with relativistic corrections subtracted. 



Hsue |62j the relativistic effects in the 3d transition metal 
series are important. Therefore, in comparing our non- 
relativistic calculations with the experiment we take such 
effects into account by subtracting the estimated values 
reported by Martin and Hay from the experimental val- 
ues. After this correction (see Table [I] for details), the 
ground state remains 3 D, with 3 F state 0.39 eV higher, 
and X S state 1.53 eV above the ground state. These val- 
ues will be referred to as "relativistically corrected (RC) 
experimental values." 

In Table|Uwe choose to utilize the Martin and Hay [6(j 
relativistic corrections as opposed to the ones computed 
by Jeng and Hsue because they include the addi- 
tional 1 S'(3(i 10 ) configuration. The results of the recent 
relativistic calculations in the RESC approximation (rel- 
ativistic scheme by eliminating small components) re- 
ported by Yanagisawa et al. |27j do not lend themselves 
to an analysis of relativistic corrections. Moreover, these 
calculations seem to be at odds with the two previous 
calculations. 

Our results, summarized in Tablc[I] show that only the 
DFT/Wachters+f calculations with B3LYP and Becke98 
hybrid exchange-correlation functionals predict a 3 D 
ground state, although B3LYP/ECP predicts the 3 D 
state only 0.01 eV above the 3 F ground state. 

It is worth mentioning that, for all our DFT calcu- 
lations, there are differences between the components of 
the 3 D state of Ni and these differences range from 4 meV 
to 37meV. We report the energy of the 3 D component 
with the lowest energy as the energy of the 3 D state. It 
is also worth mentioning that the B3LYP/ECP calcula- 
tions fail to converge for the spin a d xy -, d yz -, d xz -, and 
d x 2 _ y 2 -hole components of the 3 D state. 

The all-electron calculations with B3LYP and Bccke98 
XC functionals also predict an ordering of the 3 D, 3 F 
and 1 S states in agreement with the experiment. 

The values of the computed energies of 3 F (relative 
to D) differ from the observed experimental values by 
0.30 eV (B3LYP) and 0.26 eV (Becke98). However, when 
compared with the relativistically corrected experimen- 
tal values, the differences drop to only — 0.06eV and 
— O.lOeV, respectively. On the other hand, the com- 
puted energies of X S (relative to 3 D) are larger than the 
observed experimental values by 0.16 eV (B3LYP) and 
0.04 eV (Becke98), and larger than the relativistically 
corrected experimental values by 0.37 eV and 0.25 eV, re- 



spectively. However, the larger errors in the 1 S is less 
important for the purpose of nickel cluster calculations. 

Hartree-Fock calculations predict 3 F ground state, 3 D 
1.44 eV higher and 1 S 5.81 eV above the ground state in 
good agreement with numerical HF calculations of Mar- 
tin and Hay |60| . but with large errors compared to the 
RC experimental values. MP2 calculations predict 1 S 
ground state, with 3 D and 3 -F states 0.27 eV and 1.41 eV 
higher, respectively. 

The unoptimized FSLYP functional is, as expected, 
the least accurate. With the Wachters+f basis it yields 
results that differ from the RC experimental values and 
B3LYP and Becke98 results by « -0.5 eV for 3 F and by 
w 0.5 eV for 1 S. 

The effective core potentials (ECP) tend to overstabi- 
lize 3 F by 0.2 - 0.5 eV and destabilize X S by 0.2 - 0.4 eV 
(relative to 3 D) with respect to the all-electron counter- 
parts. Thus, all our DFT/ECP calculations predict 3 F 
ground state. However, the B3LYP/ECP calculations 
yield 3 D only 0.01 eV above the 3 F ground state, which 
can be considered acceptable error for the dissociation 
energy of nickel dimer which is of order of 2eV, given 
the savings of using ECPs. 

B. Nickel Dimer 

The determination of the ground state of N12 has been 
debated over the last few decades. According to the re- 
cent results [63l l65| . the most plausible candidates are 
spin-orbit coupled states of £1 = 0+ (a mixture of 3 E~ 
and X S+) and f2 = 0~ (a mixture of 3 S+ and 1 S„ ). 

The bond lengths (d e ), dissociation energies (-D e )|z3 
and vibrational frequencies (u> e ) for the ground state of 
N12 from different calculations are reported in Table [D] 
along with experimental values and results from other 
theoretical studies. The results in Table [H] are listed in 
the order of decreasing average absolute relative devia- 
tions (AARD) from experimental values of bond length 
(d e ), dissociation energy (D e ) and vibrational frequency 
(w e ). 

Please note that our calculations are non-relativistic 
and do not include spin-orbit coupling, and spin-orbit 
deperturbed values of molecular properties of interest for 
N12 are not available in the literature. To account for 
that, we have subtracted the CASPT2 relativistic correc- 
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TABLE II: Ground state of N12 - comparison between computations and experiment. The reported singlet states from our 
calculations are projected. d e - bond length (A), D e - dissociation energy, relative to ground state Ni atoms (without zero- 
point correction, eV), ui e - vibrational frequency (cm -1 ). The relative deviations from the experimental values are given in 
parentheses, and the average (AARD) and maximum (MARD) absolute relative deviations from experimental values of d e , D e 
and LJ e are listed under AARD and MARD columns, respectively. 



Method 


State 


de 






D e 




We 


AARD 


MARD 


FSLYP/ECP 


1 (d A d B ) 


2.236 


(1.5) 


1.325 


(-28.4) 


283.0 


(14.9) 


14.9 


28.4 


FSLYP/AE 


r (d A d B ) 


2.260 


(2.5) 


1.664 


(-10.1) 


271.1 


(10.1) 


7.6 


10.1 


Becke98/AE 


IfjA jB \ 
X^xy^xy ) 


2.296 


(4.2) 


2.071 


(11.9) 


256.8 


(4.3) 


6.8 


11.9 


CASPT2" 


ly+ 1-p 


2.281 


(3.5) 


1.89 


(2.2) 


281.0 


(14.1) 


6.6 


14.1 


CASSCF/IC-ACPF'' 


1 9 


2.291 


(3.9) 


1.691 


(-8.6) 


253.0 


(2.8) 


5.1 


8.6 


Becke98/ECP 


(d x 2_ y 2d x 2_ y 2) 


2.278 


(3.4) 


1.792 


(-3.1) 


265.1 


(7.7) 


4.7 


7.7 


B3LYP/ECP 


(d x 2_ y 2d xy ) 


2.271 


(3.0) 


1.851 


(0.1) 


269.3 


(9.4) 


4.2 


9.4 


B3LYP/AE 


(d x 2_ y 2d x 2_ y 2) 


2.291 


(3.9) 


1.835 


(-0.8) 


258.9 


(5.2) 


3.3 


5.2 


Exp. c 


0+/0" 


2.204 




1.85 




246.2 








a We report here the values from Table VIII of Ref. l63t last column (+3s3p for d e and aj e , 


and BSSE for D e ), from which we 


subtract the 



estimated relativistic corrections (RC) and, for D e only, the estimated spin-orbit coupling contributions (SO). From the same table we 
estimate the relativistic corrections to d e , D e and ui e as the difference between the values in the +RC column and ones in the CASSCF 
column, and the spin-orbit coupling contribution to D e as the difference between the value in the +SO column and the one in the +3s3p 
column. We als o subtract these RC and SO contributions from the experimental values. 
6 from Ref. 153 

c Experimental values from which we subtract the CASPT2 estimates (see footnote a) from Ref. [^3 for relati vist ic contributions (RC) for 
d e , D e and ui e , and spin-orbit contributions (SO) for D e . The experimental val ue o f d e is 2.1545 ± 0.0004 A l65l from which we subtract 
the CASPT2 RC of -0.05 A. The experimental value of D is 2.042 ± 0.002 eV from which we subtract the CASPT2 RC of 0.07 eV 
and CASPT2 SO of 0.14eV; we report D e = Dq + ^hui e . The experimental value of u> e is 259.2 ± 3.0cm — 1 (§6] from which we subtract 
the CASPT2 RC of 13cm" 1 . An earlier work [gj reported 280 ± 20cm -1 . 



tions (RC) to d e , D e and uj e , and spin-orbit contributions 
(SO) to D e from the experimental values. We estimate 
the relativistic and spin-orbit coupling corrections from 
Ref. Ib^L Please see footnote a of Table ITT1 for details. 

The reported singlet states from our calculations are 
spin-projected by the approximate method described in 
Section [H] (Computational details). 

For the results from FSLYP/ECP and Becke98/ECP 
computations, the splitting between the (d^ y d^ y ) and 
(d^2_ y 2d^ 2 _ y 2) states, both for triplet and for mixed 
S z = 0, is larger (8meV for FSLYP/ECP and 4meV 
for Becke98/ECP) than the accuracy of the DFT cal- 
culations (better than 0.1 meV). Thus, our approximate 
spin and symmetry projections are questionable for these 
particular calculations. However, since we observed even 
larger differences between the components of the 3 D state 
of Ni (up to 0.03 eV), we chose not to investigate this 
matter any further. In these cases, the reported values 
are those of the component with the lowest total energy 
(largest dissociation energy). 

Almost all of our calculations: Becke98 and FS- 
LYP both AE and ECP and B3LYP/AE predict 
1 (d^ 2 _ y2 d^ 2 _ y2 ) or 1 {d^ y d^ y ) ground state, which is 
a mixture of : E+ and l T g . B3LYP/ECP predicts a 
3 (d^2_ y 2d^ y ) (mixture of 3 £~ and 3 r„) ground state 

virtually degenerate with the l {d^_ y2 dB 2 _ y2 )/ l (d£ y d% y ) 
state, which is only 1 meV higher in energy than 
3 (d^2_ y 2d^ y ) ground state. 

Among the high-level wavefunction methods, CASPT2 
without spin-orbit coupling predicts 1 E^J" ground state 



degenerate with l T g , and CASSCF/IC-ACPF predicts 
1 r s ground state. Our DFT all-electron calculations can 
be consistent with either one of the wavefunction meth- 
ods. The experimental results are consistent with any of 
the predictions of our DFT calculations and CASPT2, 
but not with the state predicted by CASSCF/IC- 
ACPF. 

The absolute relative deviations from the experimental 
values of computed bond lengths (d e ), dissociation ener- 
gies (D e ) and vibrational frequencies (u) e ) for N12 are 
plotted in Fig. ^ arranged from left to right in order of 
decreasing total absolute relative deviation (TARD) - 
the sum of absolute relative deviations from the experi- 
mental values of the computed d e , D e and u> e . 

From Fig. as well as from Table [H] it is apparent 
that overall, for Ni2 the all-electron DFT calculations 
with B3LYP functional give the best agreement with ex- 
periment (9.9% TARD). B3LYP/ECP (12.5% TARD) 
and Becke98/ECP (14.2% TARD) follow with an over- 
all performance just a little better than CASSCF/IC- 
ACPF (15.3% TARD). Becke98/AE (20.4% TARD) 
and FSLYP/AE (22.7% TARD) are next among our 
DFT calculations, performing just a few percent worse 
than CASPT2 (19.8% TARD). With 44.8% TARD, the 
FSLYP/ECP calculation gives the largest disagreement 
with experiment and the other methods. 

The relative deviations from the experimental values of 
the computed bond length (d e ), dissociation energy (D e ), 
asymptotic dissociation energy (D®, vide infra) and vi- 
brational frequency (oj e ) of N12 are plotted in Fig. [21 for 
comparison. The values are arranged in order of increas- 
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FIG. 1: The absolute relative deviations from experiment of 
computed dissociation energy, bond length, and vibrational 
frequency for N12. The results are arranged from left to 
right in order of decreasing total absolute relative deviation 
(TARD) — the sum of absolute relative deviations from the 
experimental values of the computed d e , D e and cu e . 



FIG. 2: The relative deviations from experimental values 
of the computed bond length (d e ), dissociation energy (D e ), 
asymptotic dissociation energy (D^, see text for definition) 
and vibrational frequency (uj e ) of N12. Only the results from 
our DFT calculations are connected by lines. CASPT2 and 
CASSCF/IC-ACPF are included for comparison. 



ing deviation in the bond length. 



1. Bond length 

It is apparent that all calculations included in Ta- 
bic [H] and Fig. |2] — both our DFT calculations and 
the CASPT2 |63| and CASSCF/IAACPF wavefunc- 
tion methods included for comparison — overestimate 
the bond length of N12. The deviations from the ex- 
perimental value of the computed bond length, Ad e = 
d com P _ d ex P range between 0.03 A (1.5%) and 0.09 A 
(4.2%). 

Among our DFT calculations, the best agreement 
with the experiment for the bond length of M2 is ob- 
tained by FSLYP/ECP with Ad e = 0.032 A (1.5%), fol- 
lowed by FSLYP/AE with Ad e = 0.056 A (2.5%) and 
B3LYP/ECP with Ad e = 0.067 A (3.0%). Becke98/ECP 
with Ad e = 0.074 A (3.4%) performs very similar to 
CASPT2, for which Ad e = 0.077 A (3.5%). Both 
B3LYP/AE and CASSCF/IC-ACPF are among the 
methods that give the largest disagreement with the 
experiment, with Ad e = 0.087 A (3.9%). Finally, 
Becke98/AE yields the worst deviation from experiment, 
Ad e = 0.092 A (4.2%). 

Among all three XC functionals, the best agreement 
with experiment for the bond length is obtained with the 
FSLYP functional, both AE and ECP. B3LYP follows 
with a bond length 0.03 A longer than the one computed 
with FSLYP. Becke98 bond length is in the worst agree- 
ment with the experiment, but only w 0.005 A longer 
than the B3LYP bond length. 



For each of the three XC functionals used, ECP calcu- 
lation predicts shorter bond length than the AE one by 
« 0.02 A, and, thus, it is in better agreement with the 
experiment. 



2. Dissociation energy 

The computed dissociation energies span a large range 
of values, from 1.33 eV for FSLYP/ECP to 2.07 eV for 
Bcckc98/AE. The deviations from the experimental value 
of the computed dissociation energy, AD e = Df mp — 
D cxp range between -0.525 eV (-28.4%) and 0.221 eV 
(11.9%). 

Among our DFT calculations, the best agreement with 
the experiment for the dissociation energy of Ni2 is ob- 
tained with the B3LYP functional. B3LYP/ECP slightly 
overestimates D e by 0.001 eV (0.1%), while B3LYP/AE 
slightly underestimates D e by 0.015 eV (0.8%). This ex- 
cellent agreement with the experiment of the B3LYP 
functional is clearly fortuitous since the errors in the 
B3LYP dissociation energies average 0.10 eV, with a 
maximum absolute deviation of 0.36 eV for the G2 set of 
molecules 68]. Becke98/ECP comes second and underes- 
timates D e by 0.058eV (3.1 %), performing only slightly 
worse than CASPT2, which overestimates D e by 0.040 eV 
(2.2%). FSLYP/AE is next and it underestimates D e 
by 0.186 eV (10.1%) similar to CASSCF/IC-ACPF, for 
which AD e = -0.159 eV (-8.6%). Becke98/AE and FS- 
LYP /ECP are the methods that give the largest disagree- 
ment with the experiment: Becke98/AE overestimates 
D e by 0.221 eV (11.9%), and FSLYP/ECP underesti- 
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mates D e by 0.525 eV (28.4%). 

The effects of ECP and XC functionals on the dissoci- 
ation energy of Ni 2 do not seem to show similar trends 
to the ones seen for the bond length. However, simi- 
lar trends can be noticed if, instead of D e , one com- 
pares the asymptotic dissociation energy, D°, which is 
the dissociation energy with respect to the 3 D atoms 
that correlate with the ground state of the nickel dimer 
(D" = D e + 2E m 3 D , where E m a D is the energy of the 
3 D state of Ni atom relative to the energy of the ground 
state). 

The agreement of computed D" with the experimental 
value is clearly better than that of D e . B3LYP/AE with 
AD" = -0.015 eV (-0.8%) and B3LYP/ECP AD" = 
0.031 eV (1.7%) give the best agreement with the exper- 
iment, similar to CASPT2, for which AD a e = 0.040 eV 
(2.2%), and FSLYP/AE AD" = 0.052 eV (2.8%). FS- 
LYP/ECP with AD" = 0.106 eV (5.7%) is a little worse 
than FSLYP/AE. Becke98 gives the largest overestima- 
tion for D a e : Becke98/AE gives AD" = 0.221 cV (11.9 %) 
and Becke98/ECP gives AD" = 0.349 eV (18.8%). 

For all three functionals, the ECP basis tends to overes- 
timate the D" compared to the AE basis. For B3LYP and 
FSLYP the effect of ECP on D" is the smallest among 
the three functionals (~ 0.05 eV), while for the Becke98 
functional the effect of ECP on D" is largest (0.13 eV) , for 
which AD" increases from 0.22 eV for AE to 0.35 eV for 
ECP. However, for Becke98/ECP AD e is only -0.06 eV 
due to cancellation of large and positive AD" and the 
large D(Ni 3 D). For FSLYP this cancellation doesn't 
happen and both FSLYP/AE and FSLYP/ECP under- 
estimate the dissociation energy by fairly large amount 
because of the large error in E( 3 D Ni). 



3. Vibrational frequency 

As can be noticed in Fig. [21 there seem to be a general 
trend for all our DFT calculations, that the error in vi- 
brational frequency decreases as the error in bond length 
increases. CASSCF/IC-ACPF is close to following the 
same trend, but CASPT2 is clearly an outlier. 

It is apparent that all calculations included in Ta- 
ble [H] and Fig. |2] — both our DFT calculations and 
CASPT2 and CASSCF/IAACPF wavefunction methods 
included for comparison — overestimate the vibrational 
frequency of N12 . The deviations from the experimental 
value of the computed harmonic vibrational frequency, 
Aw e = w™ m P - co c e x P range between 10.6 cm" 1 (4.3%) 
and 36.8 cm" 1 (14.9%) among our DFT results. 

Becke98/AE with Aw e = 10.6 cm" 1 (4.3%) and 
B3LYP/AE Auj e = 12.7cm" 1 (5.2%) give the best 
agreement with the experiment among our DFT results, 
slightly worse than CASSCF/IC-ACPF, for which Auj e = 
6.8cm" 1 (2.8%). Becke98/ECP follows, overestimating 
w e by 18.9cm" 1 (7.7%). B3LYP/ECP and FSLYP/AE 
perform similarly with Au e = 23.1cm" 1 (9.4%) and 
Auj e = 24.9cm" 1 (10.1%), respectively. FSLYP/ECP 



with Auj e = 36.8 cm -1 (14.9%) gives the worst agree- 
ment with experiment, similar to CASPT2, which over- 
estimates uj £ by 34.8cm" 1 (14.1 %). 

4- Summary of the results for dfdf -holes states of Ni® 

All calculations predict dfdf -holes states to have the 
lowest energy both for singlet and for triplet spin multi- 
plicities. The bond lengths of optimized geometries, dis- 
sociation energies and vibrational frequencies for these 
states calculated with the described DFT methods are 
tabulated in Table ILTT1 for comparison. 

The first observation is that the 3 (d\ , 2 d% , 2 ) 

x y x y 

and/or 3 (d^ y d^ y ) are the highest-lying states, for all cal- 
culations, and that the spin projection changes the or- 
dering of the singlet states for all three all-electron cal- 
culations. For these calculations, the lowest energy is 
obtained for the unprojected singlet 1 ' 3 {d^ 2 _ y2 d^ y ) state, 
and upon projection, the degenerate 1 (d^ 2 _ y2 d^ 2 _ y2 ) and 

1 {dxydxy) become the ground state. 

For the B3LYP/ECP calculation spin projection does 
not change the 3 {d^ 2 _ y2 d^ y ) ground state, although it 
makes the 3 {d^ 2 __ y2 d^ y ) ground state nearly degener- 
ate with the degenerate 1 {d^ 2 _ y2 d^ 2 _ y2 ) and 1 {d^ y d^ y ). 
However, the 3 (d^ 2 _ y2 d^ y ) ground state is only 0.001 eV 
lower in energy than the degenerate 1 (d^ 2 _ y2 d^ 2 _ y2 ) and 
1 {dxy d xy)- For Becke98/ECP the unprojected ground 
state is 1 ' 3 (d^ 2 ^ y2 d^ y ) degenerate with 3 (d£ 2 _ y2 d^ y ) , 
and upon spin projection, 1 (d\_ 2 d^ 2 _, 2 ) becomes the 

v x y x y 

ground state, with a dissociation energy larger than the 
one of 1 {d^ y d^ y ) by 0.005 eV. For FSLYP/ECP the 
1 '' 3 (d^ y d^ y ) unprojected ground state does not change 
upon spin projection, but the difference between the D e 
of 1 {d^ 2 _ 2 d^ 2 _ 2 ) and that of 1 (d^ y d^ y ) is the largest 
among all calculations: 0.008 eV, and is larger than the 
numerical accuracy of the DFT calculations (better than 
0.1 meV). 

It is also worth noting that for all calculations the av- 
erage D e of singlet states is larger than the one of the 
triplet states. However, the difference between the sin- 
glet and the triplet is very small for Becke98/ECP and 
B3LYP/ECP (0.006 eV and 0.003 eV, respectively). For 
the other calculations, the difference is somewhat larger, 
around 0.015 eV. 

However, it is important to note from Table II I II that 
for each combination of exchange-correlation functional 
and basis set used, all <5(5-holes states are in a very narrow 
energy range: ~ 20meV for all all-electron calculations, 
26meV for FSLYP/ECP, 13meV for Becke98/ECP and 
only 7meV for B3LYP/ECP. 

Since, as shown above, the ordering of states can 
change upon spin-projection, if possible to perform, spin- 
projection is desirable. However, we want to emphasize 
that the differences between the lowest broken-symmetry 
singlet states and the projected singlet ground states, 
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TABLE III: DFT results for Ni 2 . d e - bond length (A), D e - 
dissociation energy, relative to ground state Ni atoms (with- 
out zero-point correction, eV), ui e - vibrational frequency 
( cm -1 ). The notation used for the states is M (h A h B ), where 
M is the multiplicity, h A and h B are the holes on Ni atoms A 
and B, respectively. The S z = 0, (S 2 ) = 1 mixed states are 
denoted by 1 ' 3 {h A h B ). 
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for the all-electron calculations and FSLYP/ECP, is less 
than 10 meV, and for some applications that difference 
may not be relevant. Nevertheless, we plan to consider 
spin-projection for larger clusters, if possible, at least for 
evaluating the errors that arise from it. 



5. Potential energy curves (PEC) 

In order to determine the ground state of M2 we did 
a full scan of the PEC for each method and for each 
unique combination of holes. All calculations predict 
<5<5-holes states to have the lowest energy, with the next 
level 50-100 meV above, a 5 for Becke98 and B3LYP cal- 
culations and irS for FSLYP calculations. 

The computations of air states with Becke98 and 
B3LYP functionals only converge to 10 _5 -10 -4 Hartree 
within 100 iterations in the 1.95-2.55 A range. Because 
the FSLYP calculations, which converge properly, pre- 
dict that these states are w 200 meV higher, the same 
value as the "not-so-converged" results for the above cal- 
culations, we have chosen not to investigate the matter 
any further. 

Since the results of the PEC scans are rather similar, 
and B3LYP is our functional of choice, in the following 
discussion of the PEC's, we focus attention principally 
on the results from B3LYP calculations. 

The B3LYP/AE and B3LYP/ECP potential energy 
curves (PEC) of singlet (S z = 0) and triplet (S z = 1) 
states of Ni2 (both unrestricted, symmetry broken) are 
shown in Fig.|21and in Fig. 21 respectively, along with the 
variation of (S 2 ) with the bond length for all possible po- 
sitions of holes in the 3d shell on both atoms, grouped 
by hole type. The first trend that can be noticed is that 
the equilibrium bond length increases as the dissociation 
energy decreases. Aside for a few states (singlet aa and 
7T7r), all states have (S 2 ) ~ 1 over a large interval, vali- 
dating the weakly interacting 3c? 9 units model for a large 
range of bond lengths. Even the singlet a a and 7T7t states 
have (S 2 ) ~ 1 in a range of about ±1 A around the equi- 
librium bond length. 

One can notice a big difference between AE and ECP 
PEC's: ECP PEC's branch around 3.5 A. There are 
two causes for branching: one, which is not related 
to functional [l^ or ECP, is the restricted-unrestricted 
crossover, while the other cause is dissociation into 3 F 
ground state of Ni atoms. These two effects overlap be- 
cause the branching is obtained by scanning the PEC 
from ps 3.5 A, increasing the bond length and using as 
initial guess the molecular orbitals from the previous cal- 
culation. Depending on the initial guess, the calcula- 
tion may end in the restricted or unrestricted solution, 
or, at large distances, the calculation may converge to 
the 3 .F + 3 F, 3 F + 3 D or 3 D + 3 D states of the Ni 
atoms. The restricted-unrestricted branching is likely 
to show up for any of the methods, but the ground- 
state branching can only appear for the methods that 
predict 3 F ground state for Ni, namely FSLYP/AE FS- 
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LYP/ECP, and Becke98/ECP along with the discussed 
B3LYP/ECP. 

One can also notice that some of the B3LYP/ECP 
PEC's have asymptotes below 0, i.e., below the en- 
ergy of the ground state of the two nickel atoms. A 
closer look reveals that the asymptotes of the d xy d xz -, 
d xy d yz - , d xz d yz - , d xz d xz -, d yz d yz -, d X yd xy - , d x 2_ y 2d xz -, 
d x 2_ y 2d yz -, d x 2_ y 2d xy -, and d ;E 2_ y 2ci : j.2_ a 2-holes states lie 
0.1 eV below the ground state of the two nickel atoms, 
the ones oid z 2d xz -, d z 2d yz -, d z 2d xy -, and d z 2d x 2 _ y 2-\\o\es 
states lie 0.04 eV below the ground state of the two nickel 
atoms, and only d^d^-holes state lies 0.03 eV above the 
ground state of the two nickel atoms. The most likely ex- 
planation for this observation is that B3LYP/ECP pre- 
dicts a lower energy for a state that is not in the space 
of states spanned by our initial guess. This issue needs 
further investigation, but since the effect is rather small 
(at most 0.05 eV/nickcl atom), we chose to investigate 
the issue in a further paper. 

The initial PEC scans are done either with broken sym- 
metry atomic initial guess (3d 9 4s 1 f f + J.J. 3d 9 4s 1 for 
singlet and 3<i 9 4s 1 f f + JT 3<i 9 4s 1 for triplet) at each 
bond length or, starting from 10 A and decreasing the 
bond length and using as initial guess the molecular or- 
bitals at the previous bond length. Either initial guess 
gives the same results, but the method using atomic ini- 
tial guess needs a few extra iterations. For larger cluster 
calculations it may be useful to save the molecular or- 
bitals at each geometry configuration and try to reuse 
them for a neighboring point calculation. 

It is apparent from Fig. that for the B3LYP/AE 
calculation the singlet dissociates to the correct 2 3 D 
atoms limit ((S 2 ) = 2), whereas the triplet dissociates 
to 3 D + 1,3 D, which is 0.14 eV above the correct limit. 
This type of error only plays an important role at large 
distances, when the molecule starts to resemble two sep- 
arated atoms, and can be correlated with (S 2 ) of the 
Kohn-Sham determinant. When (S 2 ) is close to the ex- 
act value, this type of error is not present. For Ni2, 
both singlet and triplet, the (S 2 ) is correct (i.e. equal 
to the theoretical value) for 2 A < d e < 3 A. At inter- 
atomic distances greater than approximately 3 A, (S 2 ) 
starts to increase, and so does the error in the energy of 
the triplet. At interatomic distances larger than approx- 
imately w 4 A, (S 2 ) for the triplet reaches a value of « 3 
and stays constant for larger distances. Similarly, the er- 
ror in the energy of the triplet approaches the asymptotic 
value of 0.14 eV. 

In larger clusters, this could be a potential issue for 
computing barriers. However, only configurations in 
which one atom is at sufficiently large distance from 
other atoms, completely or partly detached (evaporated) 
from the cluster, and in the 1,3 D state, would en- 
counter the above described problem. Moreover, the er- 
ror (< 0.14 eV) could be important if the height of the 
barrier were small. But the evaporation energy of an 
atom from the cluster is likely to be of the same or- 
der of magnitude as the dissociation energy of the dimer 




FIG. 5: The absolute relative deviations from experiment of 
computed dissociation energy (D e ), bond length (de), vibra- 
tional frequency (u e ), and dipole moment (fj,) for NiH. The 
results are arranged from left to right in order of decreasing 
total absolute relative deviation (TARD) — the sum of ab- 
solute relative deviations from the experimental values of the 
computed d e , D e , ui e and fi. 

(w 1.5 eV), and the height of the barrier would be over- 
estimated by w 10 %. Consequently, this error should be 
unimportant for large clusters. 



C. Nickel Hydride 

The bond lengths (d e ), dissociation energies (D e ), vi- 
brational frequencies (uj e ) and dipole moment (n) for 
ground states of NiH from different calculations are re- 
ported in Table IIVI along with experimental values and 
results from other theoretical studies, listed in the or- 
der of decreasing average of absolute relative deviations 
(AARD) from the experimental values of the computed 
d e , D e , oj e and fj,. 

The experimental values reported in Table HVI are the 
deperturbed values of d e and ui e of Gray et al. I70J. cited in 
Ref. l&j. th e recommended value of D e from Re f.l7ll (cited 
in Ref . I63T) and \i from Ref . [f^ (cited in Ref . 1631) , from 
which we subtract the CASPT2 relativistic corrections 
(RC) to d e , D^i and fi from Ref. EH and the MRCI RC to 
w e from Ref. l69l (see footnote a of Table llVl for details). 

The absolute relative deviations from the experimental 
values of the computed bond lengths (d e ), dissociation 
energies (D e ), vibrational frequencies (u e ) and dipole 
moment (/i) of the ground state of NiH are plotted in 
Fig. [S] for comparison. They are arranged from left to 
right in order of decreasing total absolute relative devia- 
tion (TARD) — the sum of absolute relative deviations 
from the experimental values of the computed d e , D e , u) e 
and [i. 
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TABLE IV: Ground state of NiH - comparison between computations and experiment. d e - bond length ( A), D e - dissociation 
energy, relative to ground state Ni atoms (without zero-point correction, eV), u e - vibrational frequency (cm -1 ) and fi - dipole 
moment (Debye). The relative errors with respect to experimental values are given in parentheses, and the average (AARD) 
and the maximum (MARD) absolute relative deviations from experimental values of d e , De, uj e and fj, are listed under the 
AARD and MARD columns, respectively. 



Method 








D e 


We 






M 


AARD 


MARD 


B3LYP/ECP 


1.454 


(-1.6) 


2.901 


(13.8) 


1937.6 


(-0.2) 


2.29 


(-12.7) 


7.1 


13.8 


Becke98/ECP 


1.456 


(-1.4) 


2.808 


(10.1) 


1927.6 


(-0.7) 


2.43 


(-7.2) 


4.8 


10.1 


B3LYP/AE 


1.474 


(-0.2) 


2.856 


(12.0) 


1940.2 


(-0.1) 


2.43 


(-7.1) 


4.8 


12.0 


FSLYP/AE 


1.470 


(-0.5) 


2.681 


(5.1) 


1943.8 


(0.1) 


2.91 


(11.2) 


4.2 


11.2 


CASPT2° 


1.463 


(-0.9) 


2.76 


(8.2) 


2022.3 


(4.2) 


2.54 


(-3.1) 


4.1 


8.2 


Becke98/AE 


1.477 


(0.0) 


2.888 


(13.3) 


1944.2 


(0.1) 


2.59 


(-1.0) 


3.6 


13.3 


FSLYP/ECP 


1.449 


(-1.9) 


2.526 


(-0.9) 


1953.4 


(0.6) 


2.74 


(4.5) 


2.0 


4.5 


Exp. 6 


1.477 




2.55 




1941.3 




2.6 


(3.8) 







"We report here the values from Table VI of Ref. \£3 for d e , D e and uj e , and from Table VII for fi [PT2F(3s3p)+RC], from which we 
subtract the estimated relativistic corrections (RC). From the same reference, we estimate the RC to d e and D e from Table V and the RC 
to n from Table VII, as the difference between the PT2F+RC values and the PT2F ones. We use the MRCI RC to ui e from Ref.^ We 
also subtract these relativistic corrections from the experimental values. 

''Experimental val ues from which we subtract the CASPT2 relativistic correc tion s (RC) to d e , D e and fi from Ref. I63t and the MRCI RC 
to uj e from Ref. l69l (see footnote a). The experimental value of d e is 1.454 A l7Ct cited in Ref . l63l . from wh ich we subtract the CASPT2 
RC of —0.023 A; the experimental value of D e is 2.70eV [recommended val ue from Ref. |7Jl cited in Ref. 1631 . from which we subtract 
the CASPT2 RC of 0.15eV; the experimental value of u e is 2001.3 cm" 1 Fol cited in Ref .163. from which we subtract the MRCI RC of 
60 cm -1 ; the experimental value of fi is 2.4±0.1Debye |?2 cited in Ref. l63l . from which we subtract the CASPT2 RC of -0.22 Debye. 



From Fig. [3 as well as from Table I1VI it is appar- 
ent that for NiH, the best overall agreement with ex- 
periment among our DFT calculations is obtained for 
FSLYP/ECP (7.9% TARD), followed by Becke98/AE 
(14.4% TARD) and FSLYP/AE (16.9% TARD) similar 
to CASPT2 (16.4 % TARD). B3LYP/AE (19.3 % TARD) 
is next, similar to Becke98/ECP (19.4% TARD), and 
B3LYP/ECP (28.3% TARD) gives the largest disagree- 
ment with experiment. 

All our DFT calculations predict 2 A (6- hole) ground 
state, in agreement with the CASPT2 calculation and 
experiment. However, it is important to note, that, for 
Becke98/ECP and FSLYP /ECP results the difference be- 
tween the d^-hole and the d 2 .2_ y 2-hole components of the 
2 A state — 2meV and 5meV, respectively — is larger 
than the error of the DFT calculations (< 0.1 meV). We 
report the energy of the component with the lowest en- 
ergy as the energy of the ground state. 

It is apparent that all calculations included in Ta- 
ble llVI and Fig. [|J] underestimate the bond length of NiH. 
Bccke98/AE with Ad e = d c ° m P - d c *P = -0.0003 A 
(—0.02 %) gives the best agreement with experiment. 
The other DFT calculations and CASPT2 give signifi- 
cantly shorter bond lengths for NiH than Becke98/AE, 
but they can still be considered in good agreement 
with the experiment, giving Ad e ranging from —0.003 A 
(-0.2%) for B3LYP/AE to -0.028 A (-1.9%) for FS- 
LYP/ECP. Among all three XC functionals, the best 
agreement with experiment for the bond length is ob- 
tained with the Becke98 functional, both AE and ECP. 
B3LYP follows with a bond length 0.003 A shorter than 
the one computed with Becke98. FSLYP bond length is 
in the worst agreement with the experiment, but only 



> 

a 

O 

> 



15 
10 
5 

-5 
-10 
-15 
-20 
-25 



— r 

<> 







bond length 
vibrational frequency 
dissociation energy 



- -<$> - dipole moment 
J I I L 



<> 



35 w 

On < 



03 



-J 



P < 
— 

in 
a, 



H 
a, 

V) 

< 



Soft 

M 
cj 

<D 

EC 



: u 
w 



ca 



q 



FIG. 6: The relative deviations of computed bond length, 
dissociation energy, vibrational frequency and dipole moment 
from experimental values for NiH. Only the results from our 
DFT calculations are connected by lines. CASPT2 values are 
included for comparison. 



ps 0.004 A longer than the B3LYP bond length. For each 
of the three XC functionals used, ECP calculation pre- 
dicts shorter bond length than the AE one by « 0.02 A, 
like in the case of Nia, but this worsens the agreement 
with the experiment, unlike in the case of M2. 

The computed dissociation energies span a large range 
of values, from 2.53 eV for FSLYP/ECP to 2.90 eV for 
B3LYP/ECP. Among all DFT computations, only FS- 
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LYP/ECP underestimates D e by 0.024 eV (0.9%), and 
gives the best agreement with the experiment. All other 
DFT computations and CASPT2 overestimate D e : FS- 
LYP/AE by 5%, CASPT2 by 8%, and Becke98/ECP, 
B3LYP/AE, Becke98/AE, and B3LYP/ECP by 10%, 
12 %, 13 %, and 14 %, respectively, giving the largest dis- 
agreement with experiment. Like in the case of N12, he 
effects of ECP and XC functionals on the dissociation 
energy of M2 do not seem to show similar trends to the 
ones seen for the bond length. It can be verified that 
trends show up upon correcting D e with the energy of 
the 3 D state of Ni, but, since for NiH there is no phys- 
ical ground for that kind of correction, we chose not to 
do it. However, it is worth noting that the errors in the 
atomic energies have such large influence on the energet- 
ics of molecules. 

The differences in the theoretical harmonic vibrational 
frequencies compared to the experimental values are less 
than 1 % for our DFT calculations, while CASPT2 has 
the largest difference from the experimental value among 
the results plotted in Fig. Eland listed in Table llVl 

For NiH the dipole moment can be expected to be a 
more sensitive measure of the quality of the method [§3 j 
and a comparison of the theoretical and experimental val- 
ues of the dipole moment listed in Table IIVI and plotted 
in Fig. H3 shows that Becke98 / AE gives the best agree- 
ment, similar to CASPT2. B3LYP/AE underestimates 
the dipole moment by 7 % and FSLYP / AE overestimates 
it by a large amount (11%). ECP have a strong ef- 
fect on (j,, lowering its value by w 0.15 D (6 %), bringing 
FSLYP/ECP in closer agreement with experiment and 
worsening the agreement for B3LYP and Becke98. It is 
worth noting that Becke98 predicts a value for \i in bet- 
ter agreement with the experiment than B3LYP. Since 
fx is a one-electron property, this may be an indication 
that Becke98 gives a more accurate ground state electron 
density. 



IV. CONCLUSION 

We have used DFT with hybrid exchange-correlation 
functionals in the broken-symmetry unrestricted for- 
malism to study the electronic structure of nickel 
dimer and nickel hydride as model systems for larger 
bare/hydrogenated nickel clusters. We have examined 
three hybrid functionals: the popular B3LYP, Becke's 
newest optimized functional Becke98, and the simple 
FSLYP functional (50% Hartree-Fock and 50% Slater 
exchange and LYP gradient-corrected correlation func- 
tional) with two basis sets: all-electron (AE) Wachters+f 
basis set and Stuttgart RSC effective core potential 
(ECP) and basis set. 

For N12 , all of our DFT calculations give bond lengths 
that are within 0.1 A (5%) from the experimental value, 
and in good agreement with the high-level wavefunction 
methods CASPT2 HI and CASSCF/IC-ACPF 64]. Only 
Becke98/AE and B3LYP/AE give harmonic vibrational 



frequencies that are within 5 % from the experimen- 
tal value, similar to CASSCF/IC-ACPF. Becke98/ECP, 
B3LYP/ECP and FSLYP / AE give uj e within 10% from 
the experimental value, similar to CASPT2, and FS- 
LYP/ECP overestimates the experimental uj e by 15%. 
The discrepancies between calculated and experimental 
values of dissociation energy span a large range, be- 
tween -28% and 12%. B3LYP/ECP, B3LYP/AE and 
Becke98/ECP give values of D e that are within less than 
5% from the experimental value, similar to CASPT2. 
FSLYP / AE and Becke98/AE give values of D e that 
are a within 12 % from experimental value, similar to 
CASSCF/IC-ACPF. FSLYP/ECP gives a value of D e 
that is smaller than the experimental value by 28%. 

For NiH, all of our DFT calculations give bond lengths 
that are within 0.03 A (2 %) from the experimental value, 
and in good agreement with CASPT2|6^|. They also 
give harmonic vibrational frequencies that are within 
less than 15 cm -1 (1 %) from the experimental value, in 
better agreement with experiment than CASPT2, which 
overestimates u> e by 4%. The discrepancies between the 
calculated and the experimental values of dissociation en- 
ergy span a large range for NiH like they do for Ni2. 
FSLYP/ECP underestimates D e by 1 %, giving the best 
agreement with the experiment. All other DFT calcula- 
tions and CASPT2 overestimate D e by amounts between 
5% and 15%. For the dipole moment the deviations 
from the experimental value span the largest range: be- 
tween -13% for B3LYP/ECP and 11 % for FSLYP/ AE. 
Underestimating it by 1%, Becke98/AE gives the best 
agreement with the experiment for the dipole moment 
of NiH, similar to CASPT2, which underestimates it by 
3%. 

We also find that for Nia, the spin-projection for 
the broken-symmetry unrestricted singlet states changes 
the ordering of the states, but the splittings are less 
than lOmeV. All our calculations predict a 55-ho\e 
ground state for M2 and 5-hole ground state for NiH. 
Upon spin-projection of the singlet state of Nia, almost 
all of our calculations: Becke98 and FSLYP both AE 
and ECP and B3LYP/AE predict 1 d fL»0 or 

1 {dxydx V ) ground state, which is a mixture of 1 S+ and 
x r g . B3LYP/ECP predicts a 3 (d£,_ ya d£ y ) (mixture of 
3 Sg and 3 T U ) ground state virtually degenerate with 
the 1 {d^_ y2 d^_ y2 )/ l {d^ y d^ y ) state, which is only 1 meV 
higher in energy than 3 (cKj 2^) ground state. The 
doublet (5-hole ground state of NiH predicted by all our 
calculations is in agreement with the experimentally pre- 
dicted 2 A ground state. For M2, all our results are con- 
sistent with the experimentally predicted ground state of 
0+ (a mixture of 1 H+ and 3 S~) or 0~ (a mixture of 
and 3 S+). 

The goal of this paper is to establish what might com- 
prise a minimally reliable method for more extensive 
nickel cluster calculations. Since none of the studied 
methods gives a good agreement with experiment for all 
computed molecular properties of M2 and NiH, we de- 
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FIG. 7: Overall performance of the studied DFT methods. 
Total values of the overall discrepancy Q are plotted, as the 
total heights of the bars, along with its components (see text 
for definition). Maximum absolute relative deviations from 
experimental values for all computed molecular properties of 
M2 (A^^ x ) and NiH (A^, 1 ^.) are also shown. 

vise an ad hoc quality indicator that we name overall 
discrepancy, Q, and we calculate it with the formula: 

Q = ^£M+^£l^-^ (4) 

i i<j 

where i runs over all 7 computed molecular properties for 
Ni 2 and NiH (d e , D e and cu e of both Ni 2 and NiH, and /x of 
NiH); ti is the relative deviation from the experimental 
value of the molecular property i; i < j stands for i,j 
running over all 21 unique pairs. 

The overall discrepancy Q is the sum of two contribu- 
tions: the average discrepancy Qa, which measures the 
overall (average) deviation of the computed molecular 
properties from the experimental values, and the consis- 
tency Q £> , which measures the consistency of the meth- 
ods both when computing different molecular properties 
of the same molecule (e.g., d e and tu e of Ni 2 ), and when 
computing molecular properties for different molecules 
(e.g., d e of Ni2 and d e of NiH). For analysis, we calcu- 
late each of the indicators Q, Qa and Qd for each of 
the molecules, by partitioning Eq.^into the components 
for Ni 2 (Q Ni2 , Q™ 2 , and Q D i2 ), the components for NiH 
(Q NlH , Q2 lH : an d QS' H )' an d the mixed components of 
Qd, Q^ l2_NlH = 55- J2i<j \ e i ~ e j\ witn i running over 
the molecular properties of Ni 2 and j running over the 
ones of NiH. 

In Fig. we plot the overall discrepancy Q along 
with its components, and the maximum absolute rela- 
tive deviations from experimental values (MARD) for all 
computed molecular properties of Ni 2 (A^ 12 ^) and NiH 
(A NiH ) 

\ max ) 

Fig. reveals that B3LYP/AE gives the lowest over- 
all discrepancy (Q — 11.2%), but followed closely by 



Becke98/AE and Becke98/ECP with a value of Q larger 
than the one of B3LYP/AE by only 0.5% and 1 %, 
respectively. They are also at the same overal qual- 
ity as CASPT2, for which Q = 12.0%. FSLYP/AE, 
with Q = 14.2% is a little worse than B3LYP/AE and 
Bcckc98/AE. It is apparent from Fig. that the use of 
ECP worsen the overall agreement with experiment. The 
largest effect of the ECP's is on the results obtained with 
the FSLYP functional, increasing the value of Q by 7.1 %. 
The effect is much smaller on B3LYP, increasing Q by 
4.6 %, and negligible on Becke98 (0.5 %). 

It can be noticed that for most of the calculations in- 
cluded in Fig. the value of Q is close to the values of 
MARD for both NiH and Ni 2 . Two methods for which 
that is not the case are worth mentioning: B3LYP/AE 
and FSLYP/ECP. Both perform significantly better for 
one of the molecules than for the other, probably by ac- 
cident. B3LYP/AE performs clearly better for Ni 2 than 
for NiH, but its MARD for NiH agrees with Q, while FS- 
LYP/ECP performs much better for NiH than for Ni 2 , 
and its MARD for Ni 2 is significantly larger than Q (by 
7.1%). Thus, FSLYP/ECP is the only method that is 
not advisable to use for bare/hydrogenated nickel clus- 
ters. However, we want to emphasize that the meth- 
ods that give the best agreement with experiment and 
CASPT2, B3LYP/AE, Becke98/AE and Becke98/ECP 
are the methods of choice. 

Our results indicate that DFT, with the B3LYP (using 
the Wachters+f all-electron basis set) and Becke98 (using 
either Wachters+f all-electron basis set or Stuttgart RSC 
effective core potential and basis set) hybrid exchange- 
correlation functionals in the broken-symmetry unre- 
stricted formalism, becomes both an efficient and reliable 
method for predicting electronic structure of our model 
Ni 2 and NiH systems, although it is far from being a 
black box method. 
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TABLE V: Charge density (CD) fitting errors for the bond lengths (d e ), dissociation energies (D e ), and harmonic vibrational 
frequencies (uj e ) of Ni 2 and NiH, and dipole moment (/x) computed with B3LYP and FSLYP functionals using "Stuttgart 
RSC ECP" ECP and basis set with "Ahlrichs Coulomb Fitting" basis. The results from the calculations using CD fitting are 
reported in the "cdfit" columns, the result from the calculations not using CD fitting are reported in "nocdfit" columns, and 
the differences between the results from the calculations using CD fitting and the results from the ones not using CD fitting 
are reported under the "cdfit err" columns, with the percent relative errors in parentheses. 





D e (eV) 


tj e (cm 1 ) 


fl( Debye) 


Mol XC nocdfit cdfit cdfit err 


nocdfit cdfit cdfit err 


nocdfit cdfit cdfit err 


nocdfit cdfit cdfit err 


NiH B3LYP 1.454 1.456 0.002 (0.1) 
NiH FSLYP 1.449 1.451 0.002 (0.1) 
Ni 2 B3LYP 2.271 2.278 0.007 (0.3) 
Ni 2 FSLYP 2.236 2.241 0.005 (0.2) 


2.901 2.784 -0.117 (-4.0) 
2.526 2.406 -0.120 (-4.8) 
1.851 1.557 -0.294 (-15.9) 
1.325 0.999 -0.326 (-24.6) 


1937.6 1998.8 61.2 (3.2) 
1953.4 2015.0 61.6 (3.2) 
269.3 255.3 -14.0 (-5.2) 
283.0 268.4 -14.6 (-5.2) 


2.29 2.21 -0.08 (-3.5) 
2.74 2.65 -0.09 (-3.3) 



TABLE VI: Averages of charge density fit errors for 
B3LYP/AE optimizations and frequency calculations for the 
12 (dfdf) singlet and triplet states of Ni 2 computed with 
and without charge density fitting. d e (mA) - bond length, E e 
(mHartree) - total energy, u> e ( cm -1 ) - vibrational frequency, 
E e (meV) - relative energies with respect to the ground state 
Ni atom, AE e (meV) - relative energies with respect to the 
lowest energy state from each type of calculation. Mean - 
mean of the differences between the computations with charge 
density fitting and those without; Stdev - standard deviation 
of the differences; Max - maximum absolute difference and 
RMS - the root-mean-square of the differences. 





d e 


E e 




E e 


AE e 




mA 


mHartree 


cm -1 


meV 


meV 


Mean 


0.30 


-0.6476 


-0.28 


-4.34 


-0.07 


Stdev 


0.02 


0.0018 


0.02 


0.05 


0.05 


Max 


0.32 


0.6497 


0.30 


4.40 


0.12 


RMS 


0.12 


0.2644 


0.12 


1.77 


0.03 



APPENDIX A: ACCURACY OF CHARGE 
DENSITY FITTING 

As stated in Section [H] we use charge density fitting 
for the calculations using the all-electron Wachters+f ba- 
sis set, for which we employ the Ahlrichs Coulomb Fit- 
ting [52u53l | basis set. For evaluating the error introduced 
by charge density fitting we perform the atomic compu- 
tations with B3LYP functional and Wachters+f basis set 
with and without charge density fitting. The charge den- 
sity fitting lowers the total energies of computed atomic 
states by 2.5-3 • 10 _4 Hartree. The errors in the relative 
energies are less severe, ranging from —3.8 to 1.5 meV. 

To be cautious, we have investigated this issue fur- 
ther by comparing results of geometry optimizations 
and frequency calculations on the 12 (dfdf) states of 
N12 (six singlet, broken symmetry, and six triplet) with 
B3LYP/AE functional both with and without charge 
density fitting. The results are summarized in Table IVT1 
Although the errors in total energies are rather large (on 
the order of a little less than 1 mHartree, as can be seen 
in column labeled £? e /mHartree in Table fVl|) . they all 
have the same sign, averaging — 0.6476±0.0018 mHartree. 



Moreover, the errors in the relative energies (with respect 
to the ground state Ni atom, labeled E e /meV in Ta- 
ble IVI|) are much smaller (« 5 meV) , and again all with 
the same sign. Finally, the relative ordering of the states 
is correct, and the root-mean-square of the relative en- 
ergies with respect to the lowest energy state from each 
calculation, labeled A.E e /meV in Table ED is 0.03 meV 
with a maximum of 0.12 meV. The maximum error due 
to charge density fitting to be expected in exploring the 
PES's of larger clusters is on the order of 2-3 meV per 
Ni atom. 

As stated in Section [HI we did not use charge density 
fitting with ECP because of the large errors that resulted 
when we tried the use of Ahlrichs Coulomb Fitting basis 
in combination with Stuttgart RSC ECP. In Table IVIwc 
report the errors in the bond lengths (rf e ), dissociation 
energies (D e ), and harmonic vibrational frequencies (w e ) 
of M2 and NiH, and dipole moment computed with 
B3LYP and FSLYP functionals using "Stuttgart RSC 
ECP" ECP and basis set with "Ahlrichs Coulomb Fit- 
ting" basis. The errors in bond lengths are negligible for 
both Ni2 and NiH, but the errors in the vibrational fre- 
quencies of both M2 and NiH, diplole moment of NiH 
and dissociation energy of NiH, of the order of 5 %, are 
significat. The error in the dissociation energy of M2 is 
large, -0.3 eV (-16%) for B3LYP and -0.33 eV (-25%) 
for FSLYP. 



APPENDIX B: ACCURACY AND 
CONVERGENCE ISSUES OF DFT 
COMPUTATIONS 

The numerical integration necessary for the evalua- 
tion of the exchange-correlation energy implemented in 
NWChem uses an Euler-MacLaurin scheme for the radial 
components (with a modified Mura-Knowles transforma- 
tion) and a Lebedev scheme for the angular components. 
Table IVlTl lists the grid details for the three levels of ac- 
curacy for the numerical integration that are used in our 
DFT calculations, labeled by the corresponding keywords 
from NWChem (medium, fine and xf ine). In the same 
table we list convergence criteria used for each level of 
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TABLE VII: Details of the integration grid for the evaluation 
of the exchange-correlation energy: the number of atomic ra- 
dial (rad.) and angular (ang.) shells for Ni and H, along with 
the corresponding convergence criteria for the DFT calcula- 
tions for each level of accuracy of the numerical integration, 
in atomic units: energy (E), density (p) and orbital gradient 
(orb. grd.). 



grid 


Ni 




H 




Accuracy 




rad. ang. 


rad. 


ang. 


E 


p orb. 


jrd. 


xf ine 


160 1454 


100 


1202 


10~ 8 


10" 7 10" 


6 


fine 


130 974 


60 


590 


io- 7 


10~ 6 10" 


-5 


medium 


112 590 


45 


434 


10~ 6 


10~ 5 10" 


-4 



accuracy of the numerical integration. 

In order to assess the errors arising from numeri- 
cal integration we have performed a series of computa- 
tions using different predefined grid schemes available in 
NWChem. First, we have performed the atomic calcula- 



tions using both xf ine and fine grids. The differences 
are of the order of total energy target accuracy of the 
fine grid (rs 1.5 • 10 ~ 7 Hartree). We have also compared 
the all-electron DFT computations using B3LYP func- 
tional with fine grid against the ones with xf ine grid 
for geometry optimization and frequency calculations for 
N12, {d^2_ y 2(i^ y ) singlet and triplet states. The differ- 
ences are on the order of 10" 4 A for equilibrium bond 
length, 2 • 10~ 6 Hartree for total equilibrium energy and 
0.2 cm" 1 for vibrational frequency. We conclude that 
the fine grid is appropriate for geometry optimization 
and vibrational frequency calculations, and have used 
it in the present work. For the potential energy curve 
(PEC) scans we use the medium grid, which gives for a 
19-point B3LYP/AE PEC scan in the range 2... 3.2 A 
of Ni 2 (d^2_ y 2d^ y ) singlet an error in energy of 16/zeV 
(maximum) and 1.4 fieV (root-mean-square) with respect 
to the fine grid computations. 
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